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Abstract 

A determinantal approximation is obtained for the permanent of a 
doubly stochastic matrix. For moderate-deviation matrix sequences, 
the asymptotic relative error is of order 0(n _1 ). 

keywords: Doubly stochastic Dirichlet distribution; Maximum-likelihood pro- 
jection; Sinkhorn projection 

1 Introduction 

A non-negative matrix of order n with unit row and column sums is called 
doubly stochastic. To each strictly positive square matrix Y there corre- 
sponds a unique re-scaled matrix A such that Yy = nAij exp(/3j + jj) where 
A is strictly positive with unit row and column sums. The doubly stochastic 
projection Y \— > A can be computed by iterative proportional scaling (Dem- 
ing and Stephan, 1940), although, in this context, it is usually called the 
Sinkhorn algorihm (Sinkhorn, 1964; Linial, Samordnitsky and Wigderson, 
1998). A non-negative matrix Y containing some zero components is said to 
be scalable if there exist real numbers jj such that = nA^ exp(a^+7j) 
where A is doubly stochastic. In that case, the pattern of zeros in Y is the 
same as the pattern of zeros in A. A non-negative matrix is scalable if and 
only if per(l") > 0. 

The permanent of a square matrix A 



per(A) =J2H A ^ 
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is the sum over permutations a: [n] — > [n] of products, -Ai,<t(i) • • • An,a(n) > 
one component taken from each row and each column. If Aij = aifijBij, 
then the permanents are related by 

n 

pei(A) = per(B) JJCa^). 

i=i 

Thus, the problem of approximating the permanent of a scalable non- negative 
matrix is reduced to the approximation of the corresponding doubly stochas- 
tic matrix, which is the subject of this paper. 

It is known that the permanent of a generic matrix is not computable in 
polynomial time as a function of n (Valiant, 1979). Since exact computation 
is not feasible for large matrices, most authors have sought to approximate 
the permanent using Monte-Carlo methods. For example, Jerrum and Sin- 
clair (1989) and Jerrum, Sinclair and Vigoda (2004) use a Markov chain 
technique, while Kou and McCullagh (2009) use an importance-sampling 
scheme for a-permanent approximation. This paper develops an entirely 
different sort of computable approximation, deterministic but asymptotic. 

The determinantal formula ([I]) is an asymptotic approximation for the 
permanent of a doubly stochastic matrix B, the conditions for which are 
ordinarily satisfied if the original matrix A has independent and identically 
distributed components. The approximation is not universally valid for all 
matrix sequences, but it is valid under moderate-deviation conditions, which 
can be checked. For example, A may be symmetric. 

2 Sinkhorn projection 

The scaling algorithm has an interpretation connected with maximum like- 
lihood for generalized linear models as follows. Let the components of Y be 
independent exponential variables with means jiij such that log//jj = Pi + 'jj 
lies in the additive subspace row+col. The maximum- likelihood projection 
Y i — y pi is such that 

>j.n, /'-, - 1) = Yj y 'i >'■-• - !) = ° 

* 3 

(McCullagh and Nelder, 1989). The residual matrix Y/ft is strictly positive 
with row and column sums equal to n, from which it follows that A = 
n~ 1 Y/ fi is doubly stochastic. The residual deviance 

Dev(A) = -2^2log(nAij) 
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is a measure of total variability for strictly positive doubly stochastic ma- 
trices, taking the value zero only for the uniform matrix Aij = 1/n. The 
exponential assumption can be replaced by any gamma distribution pro- 
vided that the gamma index v is held constant. The re-scaled matrix A 
thus generated is said to have the doubly stochastic Dirichlet distribution 
DSD n (z/) with parameter v. The smaller the value of the more extreme the 
components of A, and the greater the deviance. The distribution DSD n (l) 
is not uniform with respect to Lebesgue measure in the sense of Chatter- 
jee, Diaconis and Sly (2010), but it is presumably close for large n. The 
maximum-likelihood estimate 0(A) is such that 

21og(z>) - 2^(0) = Dev(,4)/n 2 , 

where V is the derivative of the log gamma function. It is apparent from 
simulations that if A ~ DSD ra (z^), 

20(A) x log(per(nA)/n!) = 1 + O p (n _1 ), 

with variability decreasing in n. This limiting product appears to be invari- 
ably less than one for non-Dirichlet matrices. 

The scaling Y i— > A can be accomplished either by iterative weighted 
least squares, as is usually done for generalized linear models, or by the 
Sinkhorn iterative proportional scaling algorithm. For present purposes, the 
latter algorithm is preferred because it is efficient and simple to implement. 
Linial, Samorodnitsky and Wigderson (1998) provide a a modification that 
guarantees convergence in polynomial time. For moderate-deviation matri- 
ces, the time taken to compute the Sinkhorn projection is usually negligible. 

For a doubly stochastic matrix, it is known that 

< log(per(n J 4)/n!) < n 

the lower bound being attained at the matrix A = J whose components are 
all equal, Jij = 1/n. The upper bound of log(n n /n\) ~ n is attained at each 
of the permutation matrices, which are the extreme points in the set DS n of 
doubly-stochastic matrices of order n. These bounds suggest that it may be 
easier to develop an approximation for the permanent of a doubly stochastic 
matrix than the permanent of an arbitrary positive matrix. 
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3 Moderate-deviation sequences 



A sequence of square matrices X± , Xi , . . . in which X n is of order n is called 
weakly bounded if the absolute pth moment 



is finite for all p > 1. 

The set DS ra of doubly stochastic matrices of order n is convex. The 
extreme points are the n! permutation matrices, and the central point J is 
the equally-weighted average of the extremes. Given two matrices A,B£ 
DS n , the row and column totals of the deviation matrix A — B are zero. It is 
helpful to define the L 2 -norm and associated metric in DS n by the standard 
Euclidean norm in the tangent space 



Thus, the central point has norm zero, and each extreme point has norm 



Each doubly stochastic matrix has one unit eigenvalue, and all others 
are less than or equal to one in absolute value. Consequently, it is helpful 
to define the operator norm in DS n as the Hilbert-Schmidt norm of the 
deviation A — J, i.e. the largest absolute eigenvalue |A max (yl — J)|. The 
spectral gap of A is the difference Gap(A) = 1 — |A max (A — J)|. 

A sequence of doubly stochastic matrices {A n } is said to be of moderate 
deviation if the re-scaled sequence X n = n(A n — J) is weakly bounded, 
and the spectral gap G&p(A n ) > C > is bounded below by a positive 
constant. If A, B are moderate-deviation sequences, the sequence AB of 
matrix products is also a moderate-deviation sequence, which implies that 
A' A is also of moderate deviation. Convex combinations are also of moderate 
deviation. Finally, for fixed u, the random sequence A n ~ DSD n (z^) is of 
moderate deviation with probability one. 

The main purpose of this note is to show that, for a moderate-deviation 
sequence 



with error decreasing as n — > oo. One consequence is that per(n^4)/n! is 
bounded as n — > oo. Even though it does not affect the order of magnitude 




d\A,B) = Y J {A j -B ij f 



Af = d\A,J) = Y,Al-l. 



logper(nA) = log(n!) — \ logdet(7 + J — A' 



A) + 0{n- 1 ) 
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Table 1: Numerical illustration of the determinantal approximation. 







p = l 






p = 2 




n 


logper(K) 


Approx 


Error* 


logper(K) 


Approx 


Error* 


8 


7.8905 


7.8895 


1.0178 


5.7541 


5.7454 


8.7494 


10 


11.7601 


11.7595 


0.6188 


9.0627 


9.0577 


5.0177 


12 


16.0163 


16.0159 


0.4319 


12.7701 


12.7667 


3.3501 


14 


20.5955 


20.5951 


0.3293 


16.8064 


16.8039 


2.4626 


16 


25.4522 


25.4520 


0.2636 


21.1237 


21.1218 


1.9283 


18 


30.5526 


30.5524 


0.2199 


25.6868 


25.6852 


1.5778 


20 


35.8700 


35.8698 


0.1880 


30.4682 


30.4668 


1.3310 


22 


41.3831 


41.3829 


0.1636 


35.4462 


35.4451 


1.1488 


24 


47.0743 


47.0742 


0.1534 


40.6032 


40.6022 


1.0134 



*Error x 10 3 



of the error, the modified approximation 

logper(nA) ~ log(n!) - \ log det(J + t 2 J - t 2 A'A) (1) 

with t 2 = n/(n — 1) is a worthwhile improvement for numerical work. 

Table 1 illustrates the permanental approximation applied to a class 
of symmetric positive definite matrices K{j = exp(— p\ I) for points 

xi,...,x n equally spaced on the interval [0,1], so the components range 
from exp(— p) to one. The first step uses the Sinkhorn algorithm to obtain 
a doubly stochastic matrix A from K, and the second step uses the deter- 
minantal approximation (fTJ). All values shown are on the log scale. Because 
of the difficulty of computing the exact permanent, the range of ra-values 
shown is rather limited. Nevertheless, it appears that the error decreases 
roughly as l/(300n) for p = 1, and l/(40n) for p = 2, in agreement with (fTJ). 
A sequence of matrices of this type with bounded p is of moderate deviation; 
a sequence of matrices with p oc n is not. It is unclear whether a sequence 
with p = log(n) is moderate or not, but the approximation error does appear 
to decrease roughly at the rate 1/n or log (n)/rt. 

For a very special class of Kronecker-product matrices, Table 2 shows 
the relative error in the approximation for substantially larger values of n. 
These matrices have two blocks of size n/2 each, with B^ = 1 within blocks 
and Bij = p between blocks. The smaller the between-blocks value, the 
smaller the spectral gap, and the greater the approximation error. For the 
associated doubly stochastic matrix A = 2B/(n(l + p)), Table 2 shows the 
exact value log(per(n^4)/re!) together with the error in the determinantal 
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Table 2: 


Determinantal approximation error for large : 


n. 






P 


= 0.1 


P = 


0.05 


P = 


5/n 


n 


Exact 


Err x n 


Exact 


Err x n 


Exact 


Err 


20 


0.6755 


1.6543 


1.0187 


2.8616 


0.3177 


0.0238 


40 


0.6274 


1.3823 


0.9616 


3.4398 


0.5366 


0.0258 


60 


0.6140 


1.2718 


0.9274 


3.1088 


0.6900 


0.0266 


80 


0.6082 


1.2299 


0.9114 


2.8589 


0.8076 


0.0270 


100 


0.6049 


1.2082 


0.9029 


2.7229 


0.9029 


0.0272 


200 


0.5987 


1.1698 


0.8882 


2.5229 


1.2155 


0.0278 


300 


0.5967 


1.1582 


0.8839 


2.4714 


1.4065 


0.0280 


400 


0.5957 


1.1526 


0.8818 


2.4473 


1.5444 


0.0280 



approximation — log det(/ + J — A' A) /2. Once again, the empirical evidence 
suggests that the error for moderate-deviation matrices with constant p is 
0(n _1 ), and for large-deviation matrices, 0(l/(np)). 

4 Justification for the approximation 

4.1 Permanent al expansion 

Let A be a matrix of order n such that nA\ = 1 + e£, where the row and 
column totals of e are zero. If the components of e are real numbers greater 
than — 1, A is doubly stochastic and e/n = A — J is the deviation matrix. In 
the calculations that follow, e is complex-valued and weakly bounded, and 
e/n has spectral norm strictly less than one. In other words, A need not be 
doubly stochastic, or even real, but the associated sequence is assumed to 
be of moderate deviation. 

The permanent expansion of nA as a polynomial of degree n has the 
form 

per(l + e) = ]T(1 + e ^ (1) )(l + ef ] ) • • • (1 + 

= n\ (l + n ave(ef) + ^- ave«(e[e*) + ^- ave»(e[e|e* ) + ■ ■ ■ ) 

fc=0 

where n^ k =n(n— 1) ■ ■ ■ (n — k + 1) is the descending factorial, and ave"(e® fc ) 
is the average of (n^ k ) 2 products e^ 1 • • • e[ fc taken from distinct rows and 
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distinct columns. More explicitly, the term of degree k is 

k\ e{€ >~ k\n#L> 11 ik 
with summation over distinct fc-tuples ■ ■ ■ , ife) and (ri, . . . , r&). 

In the variational calculations that follow, each term such as the re- 
stricted sum (e® k ) or the unrestricted sum J2(e® k ) is assigned a nominal 
order of magnitude as a function of n. For accounting purposes, the com- 
ponents of e are of order 0(1), and each product e[e|eji. is also regarded as 
being of order one, i.e. bounded, at least in a probabilistic sense, as n — > oo. 
A sum such as ^ "(e®* 1 ), for fixed k, is of order 0({n^ k ) 2 ) = 0(n 2fc ), so the 
L -norm 

J2(4) 2 is of order 0(n 2 ). These sums could be of smaller order 
under suitable circumstances. For example, the difference (Y — ^")e® fe is a 
sum of n 2k — (n^ k ) 2 terms, and is therefore of order 0(n 2k ~ 2 ). Thus, if the 
sum of the components of e is zero, Yl e ® k = implies that the restricted 
sum e ® k is °f order 0(n 2k ~ 2 ). Likewise, if the components of e were 
independent random variables of zero mean, Y e ir = 0(n), and a similar 
conclusion holds for the restricted sum. The variational calculations that 
follow are not based on assumptions of statistical independence of compo- 
nents, but on the arithmetic implications of zero-sum constraints on rows 
and columns. For example, A may be symmetric. 

The first goal is to show that the zero-sum restriction on the rows and 
columns of e implies that Y^( e ® k ) is of order 0(n k ) rather than 0(n 2k ). In 
other words, ave"(e® fe ) is of order 0{n~ k ). In fact, the even-degree terms 
in the permanent expansion are 0(1), while the odd terms are 0{n~ l ). 
These conclusions do not hold for all doubly stochastic matrices. For ex- 
ample, if A = (pl n + J ra )/(1 + p) for some fixed p > 0, the L 2 -norm 
1 1 -A 1 1 = \/n — lp/ (1 + p) and other scalars of a similar type are not bounded. 

4.2 Zero-sum restriction 

In the term of degree three, the restricted sum is 

= ^ S 2(eD 3 + 2(6D 3 
= 4£«) 3 - 

In the first line, the restricted sum over k ^ i,j is — e«— €j, while the restricted 
sum over t of — e* is e\ + e|. Proceeding in this way by restricted summation 
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over each non-repeated index, we arrive at the following expressions for the 
restricted tensorial sums of degree two, three and four: 



£> 2 = £ =tr( e 'e) 

E e® 3 — a ri,n,n 
t Z^i h,h,h ' 

+ 3£ tt ^^ 2 + 6 £ ae 

- 6b l^ e h,h,h,h ~ i8 Z^ e u,u,n,n ~ i8 Z^ e 
, o ri,ri,r 2 ,r 2 , r. ST~^ n,r 2 ,ri,' 



11,11,12,12 

" ri,ri,ri,n 

«1,U,«2,«2 

»*2 
«2 ■ 



Without the zero-sum constraint on the rows and columns, e ® k is a sum 
over (n^ k ) 2 e-products, and thus of order 0(n 2k ). In the reduced form, each 
distinct value occurs in duplicate at least, so there are at most k/2 distinct 
values for the row indices and k/2 for the column indices. We observe 
that £ a e® 2 = tr(e'e) is 0(n 2 ), while ^ e® 3 = 4£(eD 3 is 0{n 2 ) rather 
than 0(n 3 ). The five terms in Y} e m are 0(n 2 ), 0(n 3 ), 0(n 3 ), 0(n 4 ) and 
0(ra 4 ) respectively. The final pair can be expressed in matrix notation as 
3tr 2 (e'e) + 6tr(e'ee'e). 

For any vector x = (x\, . . . , xu) in let t(x) be the associated partition 
of [k] = {l,...,k}, i.e. r(x)(r,s) = 1 if x r = x s and zero otherwise. The 
restricted sum ^ e® k can be expressed in reduced form either as a restricted 
sum or an unrestricted sum, the only difference arising in the coefficients 
as shown above for k = 4. Generally speaking, unrestricted sums are more 
convenient for computation. For general k, the restricted sums are as follows: 

xV=zy(P)m«w £ < 

P,& V.T(l)=(J 

r:r(r)=p 

= ^m{p)m{a) ^ e[, 

p,a i:r(i)>a 
r:r(r)>p 

where the outer sum extends over ordered pairs p, a of partitions of the set 
[k], and the inner sum over the row and column indices, i = . . . , i^) and 
r = (r±, . . . , rfc), which are held constant in each block. 

The coefficients m, rr$ are multiplicative functions of the partition 

™Hp) = " !)' m (P) = ]\(-V* b ~\#b - 1)!, 
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and ms(p) = m{p) = if p has a singleton block. Here, and elsewhere, 
If. denotes the maximal partition of [k] with one block, j^p is the number 
of blocks, and #6 is the number of elements of block b G p. Although 
m and are related by Mobius inversion, these expressions are not to 
be confused with the Mobius function for the partition lattice: M(p, 1) = 
(— i)# p ~ 1 (#p— 1)!, which is a function of the number of blocks independently 
of their sizes. 

For simplicity of notation in what follows, we write e p a for the sum, 

e P = e r 

i:r(i)>a 
r:r(r)>p 

in which i and r are constant within blocks, but the values for distinct blocks 
may be equal. Putting these expressions together with a scalar < t < 1, 
we find 

per(l + te) A t k ^ 

^ = E ™(p)m(ay a (2) 

k=0 cr,peV k 

where Vk is the set of partitions of [k\. 

The polynomial (|2|) is the exponential generator Y2k=o t k Ck/k\ for the 
moment numbers 

Ck = ^k E rn(p)m(a)e p a 

taking Qk = for k > n. The coefficients in the log series are the cumulant 
numbers 

c^Ec- 1 )^" 1 ^- 1 )']!^ 



= E m(p)m(aK E (-^(^ - 1)! J] ^ 

p,cr£Pk r>pVo- 6Gr 

= E ™(/'W^ xA n(/ ,Vff ) 

for fc = 1, Note that /?, <r in line 2 are the restrictions to the blocks of t 

of the partitions p, a in line 3, so r > p, a. Evidently, 

A^)=E(- i ) # ^ i (#-- i ) ! n^ 

tj>p 6GiT 
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is the generalized cumulant associated with the reciprocal descending facto- 
rial series. For example, 

r^ 4 A„(12|34) = 1 - 
n+ 5 A n (123|45) = 1 - 



n l4 


4n-6 








6(n - 2) 






3n 46 


2n l6 



n^A (12134156) - 1 - + - 8(n - 3)(7n - 10) 

n A„(12|34|5b) - 1 + ( ^ 2)3 - ( ^ 2)2 

Thus, the complete formal expansion for the log permanent is 

log(per(l + te)/n!) = - ^ m(p)m(a)e p a A n (p V a) 

k=l p,<r£Vk 
oo 

= EfciE A ™( T ) E rn{p)m{a)4 

k=l ' r£P k p\a=T 

in which A n (lfc) = l/n^ k for k < n, and zero otherwise. 

4.3 Moderate-deviation asymptotic expansion 

For a fixed pair of partitions p,a £ Vk the key scalar e„ is the sum over 
indices r(i) > a and r(r) > p of e[. Under the rules for determining the 
order of magintude in n, e p is of order 0(ra# f7+ # p ), and 

A n (p V <r) e£ = 0(n #<T+#p - fc - # ^ V(T ) +1 ). 

Since neither partition has singleton blocks, 2#a < k and 2#p < fc. Thus 
the order is 0(1) only if k is even, each partition is binary with all blocks 
of size two, and the least upper bound is the full set Otherwise, if the 
least upper bound is less than or if any block of either partition is of size 
three or more, the term is 0{n~ v ) or smaller. 

The leading term of maximal order in the expansion of the log permanent 

is 

k=l y ' p,a~2 k 



p\/a=l 2k 

\x((f'f 

2k 



n l 2 f 2k 

EW tr « e ' £ ) fe )/ n ' 2fc 
k=i 
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oo 2k 

= E^ tr (( e ' £ /^ 2 ) fc ) + 0(n- 1 ) 
k=i 

= -\ log det(J„ -t 2 e'e/n^) + 0(n~ 1 ) 
>4 0) (1) = -^logdet(/ n + J - A'A) + O^- 1 ). 

The symbol p, a ~ 2 fc denotes two partitions of [2fe] having /c blocks of size 
two, and there are (2k — 1)1 pairs whose least upper bound is the full set 
[2k]. The series is convergent for \t\ < l/|A max (A — J)|, so hn\l) is finite by 
the spectral gap assumption. 

In the expansion of the first-order correction h^(t), the terms may be 
grouped by degree in t. The following are the terms of degree six or less, 
expressed so far as possible using matrix operations in which e r and (e'e) r 
are the component- wise rth powers of e and (e'e) respectively. 

<•*-* JLv> = ^y> n 3 

degree 4: ^(e' 2 e 2 + e 2 e 2 ) 

2 \ 

degree 5: — - trfee'ee') H f / £2^62 

degree 6: ^ X^ e ' e ) 3 + ( ee ') 3 ) + 5Z( e 2 ee ' e 2 + e 2 e'ee 2 ) 

-^E( e ^) 2 + 4(«')2) 

5 Random matrices 

In order to test the adequacy of the determinantal approximation, 5000 ran- 
dom matrices A of order 10-25 were generated, and the exact values y(A) = 
log (per (nA)/n\) were computed using the Ryser algorithm. The matrices of 
even order were generated from the distribution A^ ~ DSD„(fj). In order 
to ensure an adequate range of permanental values, the parameters i/j were 
generated randomly and independently from the Gamma distribution with 
mean one on four degrees of freedom, i.e. xf/8. On the matrices of odd order 
was superimposed a multiplicative random block pattern (77 + rj'B(r,s)) in 
which 77, 7/ are uniform (0, 1) random scalars, and B is a random partition 
of [n], the Chinese restaurant process with parameter 1. This block pattern 
was superimposed prior to Sinkhorn projection in order to vary the eigen- 
value pattern, which depends on the coefficients 77 and on the block sizes. 
Matrices of this type do not satisfy the moderate-deviation condition. 



11 




s 2! „ . L 22 




~i 1 1 r 



:J§§! 



m 



i 




10 12 14 16 18 20 22 24 



10 12 14 16 18 20 22 24 



Figure 1. Top row: exact permanent of 500 random matrices plotted against 
two approximations, both on the logdog scale. Middle row: standardized logdog 
residuals plotted against the two approximations. Third row: residuals plotted 
against n, at twice the scale in the right panel. 
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On the log scale, the range of observed permanental values was 0.11 < 
y(A) < 7.56, or 0.042 < y/log(n) < 2.57, the largest value occurring for a 
matrix of order 19. The target range y(A) < 2 was exceeded by 2.6% of 
the matrices generated, and the moderate-deviation threshold y > log(n) 
was exceeded by 43 matrices comprising roughly 0.9% of the simulations. 
Most of these exceedances occurred for matrices of odd order having a pro- 
nounced block structure. The range of L 2 -values was 0.22 < ||A|| 2 < 7.32, 
with only 2% in excess of 3.0. Although there exist matrices A £ DS n such 
that the ratio y(^4)/||A|| 2 is arbitrarily large or arbitrarily small, the range 
of simulated values was only (0.53,1.15). This is a reminder that the be- 
haviour of the permanent in a bounded L 2 -ball is not a good indicator of its 
behaviour near the corners. 

Two approximations were computed, the first-order determinantal ap- 
proximation xo(A) as in (fT]), and a second-order modification x%(A) using 
the additional terms in section 4.3. The scatterplots of loglogy(^4) against 
log log xq (A) and loglogxi(j4) are shown in the top two panels of Fig. 1. To 
reduce clutter in the top two rows, only a 10% sample is shown. 

Since the relative error of the approximation x(A) increases with its 
magnitude, the standardized residuals are most naturally defined on the 
log-log scale 

Reside x) = ^oMA)/x(A)) _ 

In the middle panels of Fig. 1, the log- log residuals are plotted against x, 
using the same scale for both plots. In the lower panels, the log-log residuals 
are plotted against n. The first plot makes it clear that the determinantal 
approximation xq(A) tends to be an under-estimate: the rate of occurrence 
of the inequality y(A) > xo(A) increases from 44% for n = 10 to over 90% 
for n > 20. In the lower left panel, the alternating pattern of residuals 
for odd and even n is a consequence of the block pattern embedded in the 
matrices of odd order. For both types of matrices, the plots suggest that 
the residual distribution is asymptotically constant as n — > oo, and that 
the correction term is appreciable for moderate n. The root-mean-squared 
residual is approximately 1.2/n for the first approximation and 0.3/n for 
the second, but the distributions are more like Cauchy than Gaussian. As 
it happens, the inequalities 

logy(A) < logx (A) - lAx (A)/n, logy {A) > logx (^) + 3.0x (A)/n 
occur in the sample with rates 1.0% each. The inequality 
| logy (A) - logxi(^)| > 0.5xi (A)/n 
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occurs at a fairly constant rate of around 1% for non-block-structured ma- 
trices. For block-structured matrices, the rate is about five times as high, 
but slowly decreasing in n over the range observed. Of course, these rates 
must depend on the distribution by which the matrices are generated, but 
it is anticipated that the rate will have a positive limit for matrices in the 
moderate-deviation region xq(A) < log(n). In other words, we should ex- 
pect log(per(n.A)/n!) to lie in the interval x\(A) exp(±xi(A) log(n)/n) with 
probability tending to one for large n if x±(A) is bounded. 
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